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Abstract 

Aims: Our aim is to investigate the resistive relax- 
ation of a magnetic loop that contains braided mag- 
netic flux but no net current or helicity. The loop 
is subject to line-tied boundary conditions. We in- 
vestigate the dynamical processes that occur during 
this relaxation, in particular the magnetic reconnec- 
tion that occurs, and discuss the nature of the final 
equilibrium. 

Methods: The three-dimensional evolution of a 
braided magnetic field is followed in a series of re- 
sistive MHD simulations. 

Results: It is found that, following an instability 
within the loop, a myriad of thin current layers 
forms, via a cascade-like process. This cascade be- 
comes more developed and continues for a longer pe- 
riod of time for higher magnetic Reynolds number. 
During the cascade, magnetic flux is reconnected 
multiple times, with the level of this 'multiple re- 
connection' positively correlated with the magnetic 
Reynolds number. Eventually the system evolves 
into a state with no more small-scale current lay- 
ers. This final state is found to approximate a non- 
linear force-free field consisting of two flux tubes 
of oppositely-signed twist embedded in a uniform 
background field. 

1 Introduction 

The braiding of magnetic loops in the solar corona, 
via convective motions at the solar surface, has long 
been suggested as a potential mechanism for heat- 



ing the corona. Parker (1972) proposed that in re- 
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sponse to arbitrary footpoint motions at the photo- 
sphere, the coronal magnetic field will ideally relax 
towards a force- free equilibrium containing tangen- 
tial discontinuities, corresponding to current sheets. 
Therefore, it was argued, as such a singular state 
is approached during the relaxation, the diffusiv- 
ity of the plasma must always become important, 
leading to magnetic reconnection and plasma heat- 
ing. Since this idea was first proposed, there have 
been many arguments for and against its validity 



(e.g. van Ballegooijen, 1985; Longcope and Strauss 



1994; Ng and Bhattacharjee, 1998). 



In Parker's original 'topological dissipation' 
model, an initially uniform magnetic field is taken 
between two perfectly conducting parallel plates. 
Then random (smooth) motions are applied on the 
perfectly conducting plates. A number of attempts 
have been made to simulate such a scenario nu- 
merically, using a variety of numerical approaches. 



Craig and Sneyd (2005) employed an ideal relax- 



ation scheme which includes a fictitious frictional 



term in the equation of motion (Craig and Sneyd 



1986). They applied various complex deformations 



at the boundaries of the domain. They found that 
regardless of the nature or extent of the deforma- 
tion, no tangential discontinuities developed in the 
relaxed state, with the current concentrations re- 



maining large-scale in all cases. Mikic et al. (1989) 



employed a slightly different approach. They began 
with a uniform field and then sequentially applied 
large-scale shear flows of random orientation on one 
boundary (while the field at the other boundary was 
held fixed). In the time between each boundary 
shear, they solved the ideal MHD equations (ne- 
glecting the pressure and advective terms in the 
equation of motion), including a large spatially- 
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uniform viscosity in order to relax the magnetic 
field towards equilibrium. They found, once again, 
that no discontinuities were formed during any of 
these relaxation processes. However, progressively 
smaller scale current structures were found to de- 
velop as the number of shear disturbances increased. 
The implication is therefore that after a sufficient 
length of time, the scales would reach those appro- 
priate for dissipation in the corona. They found, 
furthermore, that the current density in the do- 
main increased exponentially in time, which is con- 
sistent with the earlier analysis of | van Ballegooijen| 



(1985), who predicted thin non-singular current lay- 



ers rather than tangential discontinuities. 



Galsgaard and Nordlund (1996) employed a fur 



ther different approach to those discussed above, 
solving the full set of resistive MHD equations with- 
out employing any artificial force (such as enhanced 
viscosity) to inhibit the plasma dynamics. They ap- 
plied a similar sequence of large-scale shearing mo- 
tions as Mikic et al. (1989), this time at both driv- 



ing boundaries. The amplitude and orientation of 
the shearing were chosen at random (from a normal 
distribution). The results were similar to those of 
Mikic et al.| ( |1989 ), in that after only a few Alfven 
crossing times, multiple small-scale current fila- 
ments were found to form in the domain. However, 
with the shearing applied on both boundaries, ex- 
ponentially growing currents were already obtained 
after the second shear motion due to the interlock- 
ing of the field lines, with the tension force setting 
up a stagnation flow. Due to the finite resistivity in 
their simulations, magnetic reconnection and Joule 
dissipation occurred in these filaments. This en- 
ergy release was found to be intermittent or 'bursty' 
when the time scale for energy input (via the bound- 
ary shearing) was long compared with the Alfven 



crossing time. The work was extended by Galsgaard 



(2002) with the inclusion of a stratified atmosphere. 
Furthermore, the generation and properties of tur- 
bulence in Parker's model have been studied in the 
framework of reduced MHD by a number of authors. 
In the initial studies, two-dimensional simulations 
were performed, with imposed 'forcing terms' tak- 
ing the place of the boundary driving (e.g. lEinaudi 



et al., 


1996 


et al. , 


1998 



|Dmitruk and G6mez[ |1997[ |GeorgouIIs 



More recently, three-dimensional sim- 
ulations of Parker's model with the reduced MHD 



2003; Rappazzo et al., 2008), and the resulting en- 
ergy spectra and heating event statistics investi- 
gated. 

In the present work, we approach the topological 
dissipation problem discussed above from a differ- 
ent angle to previous studies. We take a braided 
magnetic field that is close to force-free as an initial 
condition, and concentrate on the details of the sub- 
sequent evolution, including current sheet forma- 
tion. In general, braided force-free fields could arise 
in the solar corona either through the emergence 
of braided flux from the interior or through ran- 
dom footpoint motions at the solar surface, accom- 
panied by some intermediate relaxation processes. 
The ideal relaxation of this field towards a force-free 



state was considered in |Wilmot- Smith et al. (2009a 



where only large-scale current features were found 
in the end state. We now follow the evolution of 
the system in a resistive MHD simulation. In a pre- 



vious paper ( Wilmot-Smith et al. (2010), hereafter 



equations have been performed (e.g. Dmitruk et al. 



referred to as Paper I) we described how during the 
early evolution an instability occurs which moves 
the system away from equilibrium, leading to the 
formation of current sheets and thus the onset of 
magnetic reconnection. In the present paper we ad- 
dress the following key questions: 

1. As the plasma seeks a new equilibrium, what is 
the nature of the resistive relaxation process? 

2. In particular what are the properties of the 
magnetic reconnection processes that take 
place? 

3. What is the nature of the final state of the re- 
laxation? 

4. What is the dependence of each of the above 
on the magnetic Reynolds number? 

In Section [2] we introduce the numerical method 
and summarise the results of our previous investi- 
gations. In Section [3] we describe qualitatively the 
evolution of the system, and in Section [4] we discuss 
the nature of the final state. In Section [5] we investi- 
gate the reconnection in the system in more detail, 
and then in Section [6] discuss the dependence on the 
plasma resistivity. Finally in Section [7] we present 
our conclusions. 
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2 Numerical method and early 
evolution 

2.1 Numerical scheme 

The numerical scheme employed in the simulations 
that follow is described briefly below (further details 
may be found in Nordlund and Galsgaard| ( |1997 ) 
and at http://www.astro.ku.dk/~kg). We solve the 
three-dimensional resistive MHD equations in the 
form 



d_ 
di 



dB 

~dt 
E 

J 

dp 
di 

(Pv) 

de 
di 



-V x E, 

- (v x B) 
V x B, 



-V • (pv) , 

-V • (pvv + t) - VP + J 
-V • (ev) - P V • v + Q visc 



(1) 
(2) 
(3) 

(4) 
< B,(5) 
+ 03) 



where B is the magnetic field, E the electric field, 
v the plasma velocity, p the density, 77 the resis- 
tivity, J the electric current density, r the viscous 
stress tensor, P the pressure, e the thermal en- 
ergy, Qvisc the viscous dissipation and Q j the Joule 
dissipation. An ideal gas is assumed, and hence 
P = (7 — 1) e = |e. The equations above have 
been non-dimensionalised by setting the magnetic 
permeability po = 1, and the gas constant equal to 
the mean molecular weight. The result is that for a 
volume in which \p\ = |B| = 1, time units are such 
that an Alfven wave would travel one space unit in 
one unit of time. 

In all simulation runs we employ an explicitly 
prescribed, spatially uniform resistivity, and thus 
we have simply Qj = r]J 2 . Viscosity is calcu- 
lated using a combined second-order and fourth- 
order method (sometimes termed 'hyper- viscosity'), 
which is capable of providing sufficient localised dis- 
sipation where necessary to handle the development 



of numerical instabilities (Nordlund and Galsgaard 
1997). The effect is to 'switch on' the viscosity 
where very short length scales develop in v, while 
maintaining a minimal amount of viscous dissipa- 
tion where the velocity field is smooth. As such, the 
total viscous dissipation in the domain is negligible 
when compared with the Joule dissipation. All sim- 
ulations are carried out on a grid with 512 3 nodes. 



The dimensionless plasma density is initialised to 
have value p — 1 throughout the domain, while the 
thermal energy is initialised with value e = 0.1. The 
result is a plasma- (3 that varies throughout the do- 
main between 0.10 and 0.14 at t = 0. 

We have checked our numerical results by re- 
peating the standard run described below (with 
= 10 
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using the Lare3d MHD code (Arber 



et al., 2001). The qualitative results described in 



sections [3H6] have been verified with this different 
scheme, demonstrating their robustness. 

2.2 Initial conditions and early time evo- 
lution 

We use as an initial condition for our resistive MHD 
simulations a magnetic field that is close to being 



force- free QWilmot-Smith et al.[ |2009a[ ) , and that is 
obtained via an ideal Lagrangian relaxation scheme 
([Craig and Sne^d[ p86| |Pontin et al.[ [2009]). This 



is explained in detail in Paper I. In what follows, 
we term the input field for the ideal relaxation 
scheme the 'pre-initial' magnetic field. The final 
state of this relaxation, which is the initial condi- 
tion for the resistive MHD simulations described 
herein, is termed the 'initial' magnetic field. The 
pre-initial and initial magnetic fields have the same 
braided magnetic topology (as the relaxation is ex- 
actly ideal). The magnetic field is modelled on the 
so-called pigtail braid, with a subset of the field lines 
being linked in a non-trivial way - see Figure Qa). 
The pre-initial structure consists of six localised re- 
gions of twisted flux in an otherwise homogeneous 
field. These twists are of equal magnitude, but with 
opposite sign in three of the regions, such that the 
net current density and net helicity in the computa- 
tional volume are both zero. (To be precise we cal- 
culate the relative helicity QBerger and Field , 1984) 
using the potential field (B p = [0, 0, 1]) as the refer- 
ence field. In our case we can choose a gauge such 
that A p = A on the boundary so that the relative 
helicity reduces to an integral over A • B only.) In 
the initial magnetic field used for our MHD sim- 
ulations, there are two large-scale, diffuse current 
structures within the domain, but no thin current 
layers - see Figure [ijb) and Wilmot-Smith et al.| 
( 2009a|b ) for a full description. 

The numerical domain has dimensions [x, z] E 
[±6, ±6, ±24]. The magnetic field is line-tied onz = 
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(a) 




y 



Figure 1: (a) Three representative magnetic field 
lines from the initial magnetic field (t = 0) which 
demonstrate the braiding present in the field, (b) 
Isosurface of current density |J| at t = 0, at 25% of 
its maximum value. 



±24, where the plasma velocity is fixed at zero for 
the whole experiment, while the x and y boundary 
conditions are periodic. Since B = [0, 0, 1] at the x- 
and ^/-boundaries, these boundaries are closed with 
respect to the magnetic field at t = 0. We find that 
the dynamics is confined within the domain (x, y E 
[—4,4]) and, with the confining background field, 
field lines do not leave the numerical box through 
the side boundaries at any stage of the simulation. 
We note finally that the simulations are run for a 
sufficient time that in the final state the resistive 
dissipation is occurring on the global length scale 
of the regions of twist that were present in the pre- 
initial field (see below). 

The evolution of the system at early stages can 
be summarised as follows, and is described in detail 
in Paper I. The system initially remains relatively 
stationary, being close to a force-free equilibrium. 
However, at approximately t = 8 an instability sets 
in, and the current intensifies in two thin layers. 
At present the nature of this instability is not fully 
understood, though we propose that it may be a 
resistive instability analogous to the internal kink 
instability. This time interval corresponds to the 
Alfven travel time between two of the imposed re- 
gions of twist in the pre-initial magnetic field. We 
go on below to investigate the continued evolution 
of the system. 

3 Long time evolution 

3.1 Qualitative behaviour 

Following the instability which gives rise to the ini- 
tial current layers, the system goes on to develop 
an increasing number of thin current layers, as il- 
lustrated in Figure [2] where isosurfaces of current 
density at 50% of the domain maximum are shown. 
It is worth noting that the two initial current lay- 
ers are stronger than any of those that form in the 
subsequent evolution (see the following section and 
Figure [4]) . The development of a highly fragmented 
current structure means that magnetic reconnection 
also occurs in a highly fragmented portion of the do- 
main. This magnetic reconnection is investigated in 
detail in Section[5j The evolution resembles qualita- 
tively the onset of a turbulent-like behaviour, with 
the current layers forming a highly disordered pat- 
tern which fills the central portion (in xy) of the 
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Figure 2: Isosurfaces of the current density |J| at t = 15,27.5,35,50,290, at 50% of the corresponding 
maximum value of |J|. For the run with 77 = 10 -3 . 



computational domain (the region in which the ini- 
tial field line mapping is non-trivial). Investigating 
the power spectra of relevant physical quantities to 
ascertain whether indeed true turbulence develops 
is beyond the scope of this paper. It is in any case 
likely to require higher numerical resolution (and 
the lower value of imposed 77 that this permits) in 
order to resolve a sufficient number of decades of 
spatial scales. Thus, in the remainder of this paper 
we use the terms turbulence and cascade only in a 
very loose sense. 

In Paper I, it was shown that the locations where 
the first pair of current sheets form are well pre- 
dicted by analysing in combination the squashing 



factor Q of the field line mapping (Titov et al. 



2002; Titov, 2007) and the integrated parallel cur- 



rent along field lines, J = J J\\dl. The initial cur- 
rent layers form in regions threaded by field lines 
with high values of both quantities (for an explana- 
tion of the importance of J see Wilmot-Smith et al. 
(2009a}). The question naturally arises: do the sec- 
ondary, tertiary and subsequent current layers form 
in response to the dynamics that are set in motion 
by the initial instability and reconnection process? 
Or are they the result of multiple instabilities which 
occur throughout the domain with different growth 
rates - with the initial current sheets that we ob- 
serve being simply those associated with the fastest 





Figure 3: Current density |J| at z — for [x,y] E 
[±3.5, ±3.5] every t = 20 units from t = to t = 340 
(top-left to bottom-right). The shading intensity is 
scaled individually to the maximum at each time, 
for clarity. For the run with 77 = 10 -3 . 



growing instabilities? The myriad of thin layers of 
high Q and J found in the domain would support 
the latter idea. The evidence suggests that neither 
of these alternatives is solely responsible, but rather 
a combination of the two. However, this is another 
aspect of the simulations that warrants further in- 
vestigation. 

In an attempt to gain some insight, one can look 
at individual reconnection events in the domain. 
One can identify the outflow from certain current 
layers with the inflow to another reconnection pro- 
cess (this being one typical picture of how a mul- 
titude of current layers may form, where one re- 
connection process drives the next and so on (e.g. 
Galsgaard and Nordlund[ |1996 ; Watson and Craig , 
2003[ [Hood et al. , 2009)). However, there are many 



other current layers for which this is clearly not the 
case. Moreover, determining the local relationship 
between reconnection processes is not straightfor- 
ward - in three dimensions reconnection need not 
necessarily occur in a region of hyperbolic field, and 
is also not always associated with a hyperbolic flow 
structure with well-defined inflow and outflow. In 
fact, in three dimensions the characteristic property 
for reconnection in the absence of null points is a 
counter-rotational plasma flow on either side (with 
respect to the direction of B) of the diffusion region 
(Hornig and Priest, 2003[ ). Thus, determining the 
dynamic interaction between different reconnection 
processes within the volume is a non-trivial prob- 
lem. 

The increase in the complexity of the current dis- 
tribution continues until around t = 80, after which 
the number of current layers begins to decrease. 
Eventually, the current concentrations again attain 
large scales (see Figures [2] and [3]) . The properties 
of the final state of the simulation are considered 
in detail in Section [U The extent and duration of 
the fragmentation of the current layers - and the 
subsequent time taken to return to a configuration 
with currents on system-size scales - is found to be 
dependent on the value taken for the resisitvity, 77. 
This is discussed in Section |6l 

3.2 Relaxation properties 

The magnetic field begins in a state that is close to 
a (non-linear) force-free equilibrium. As mentioned 
above (and described in more detail in Paper I), an 
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Figure 4: From top to bottom: Time evolution of 
the spatial maximum of |J|; the Joule dissipation 
77J 2 per unit volume, the magnetic energy per unit 
volume, and the thermal energy per unit volume, 
for the run with 77 = 10 -3 . 



instability is then triggered forcing the system away 
from equilibrium. The subsequent resistive evolu- 
tion can be described as a relaxation process, which 
is mediated by magnetic reconnection (see Section 
[5]) . This relaxation process is demonstrated in Fig- 
ure [4] where the evolution of various quantities with 
time is shown. The top panel shows a time series 
of the peak current in the domain. After an initial 
sharp drop (that is likely due largely to the change 
from one discretisation to another), this soon peaks 
(at t = 12, see Paper I) before gradually declining 
towards the end of the simulation run. The current 
settles to a constant value after t « 200. After this 
point the configuration is approximately stationary, 
although the simulation has been run to t = 350 to 
be sure we are in a state whose topology is stable. 
The Joule dissipation per unit volume shows a sim- 
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Figure 5: Total per unit volume of the kinetic en- 
ergy plus the magnetic energy in excess of that of 
the background field, plotted on a log-log scale, for 
the run with 77 — 10 -3 . 



ilar pattern of evolution to the peak current. The 
free magnetic energy in the initial state is relatively 
modest, with only 0.5% magnetic energy per unit 
volume above that of the potential (homogenous) 
field. As time progresses, this magnetic energy de- 
cays, being converted into kinetic and thermal en- 
ergy (see the lower panel of the figure). We note 
that in the final state, the magnetic energy is still 
around 0.2% above that of the homogeneous field. 
In total over 60% of the free energy in the initial 
configuration is converted to other forms. The na- 
ture of the equilibrium approached is discussed in 
detail in the following section. 

As noted above the relaxation process has a tur- 
bulent appearance. However, due to the moderate 
Reynolds numbers that we are able to use with the 
present resolution, the system is unlikely to be in a 
state of fully developed turbulence, though we an- 
ticipate that such a state would be reached were we 
able to perform simulations with lower values of the 
resistivity. In terms of the categorisation described 
by |Biskamp| (2003), the system would then be de- 
scribed as undergoing decaying turbulence, since 
the additional energy in the system is included as 
excess magnetic energy in the system at large scales 
at t = (in contrast with the forced turbulence sim- 
ulations in which a mechanical stressing is continu- 



decay exponentially in time during decaying tur- 
bulence (Biskamp, 2003). In Figure [5] this quan- 
tity is plotted versus time on a log-log scale, where 
we have subtracted the constant contribution to 
the magnetic energy arising from the uniform back- 
ground field, for clarity. The dashed line shows a 
least-squares fit for the data points between t = 20 
and t = 180 when the current distribution contains 
small scales. It is clear that we do indeed see a 
power law dependence during this time, which is 
suggestive that we are at least approaching a regime 
of decaying turbulence. It is also apparent from the 
plot that after t « 180 (after which physical quanti- 
ties vary on system-size scales), this power law is no 
longer followed, with a steepening of the curve sug- 
gesting a possible transition towards an exponential 
decay. 

While the simulation has been run to t = 350, it 
is found that after t « 200 the topology of the mag- 
netic field changes only due to diffusion in the two 
tubular current structures which have system-size 
dimensions both along B and perpendicular to B 
(see Figures [2] and [3j and Figure [9] bottom panel). 
Owing to the enormous diffusive timescale in the 
corona, this diffusive evolution is not relevant for 
the coronal field. We therefore consider the state 
reached at the end of our simulations to charac- 
terise the physically relevant final state. We there- 
fore take t = 290 as the final state of the simulation 
for our investigations into its nature that follow. It 
should be noted that this final state is not an exact 
force-free equilibrium, due to the finite plasma pres- 
sure, as well the residual Alfvenic oscillation of the 
structure about its equilibrium (our use of a hyper 
viscosity means that these oscillations are damped 
only very slowly). However, it is close to force- free, 
in that in terms of the dimensionless measure of the 
proximity to a force- free field, eZ, given in Wil mot-] 
Smith et al. (2009a), we have el^6x 10 -3 (where 



ally applied, e.g 


Einaudi et al. 


(1996); 


Dmitruk and 


Gomes 


1 (1997); 


Georgoulis et al. 


(1998 


); 


Dmitruk 


et al. 


(2003); Rappazzo et al. ( 


2008|)). 


The sum 



of the magnetic and kinetic energies is expected to 



el E [0, 1] and el = for a perfectly force- free field). 

4 Nature of the final state 

4.1 Properties of the relaxation 

In order to investigate the nature of the final state, 
we take a grid of starting points on the lower bound- 
ary z — — 24 (taking 160 2 points over the region 
[—3, 3] 2 ) and trace field lines through the volume to 
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Figure 6: The mean value of a* along field lines 
shown on a section of the lower boundary (z = —24) 
of the domain in the initial (top, t = 0) and final 
(bottom, t = 290) states showing a smoothing of a* 
during the resistive relaxation with r\ = 10 -3 . 



the upper boundary. For each field line, we calcu- 
late a*, the mean value along the field line of the 
quantity a* = (J ■ B)/(B • B). Note that for a 
perfectly force-free field this simply reduces to the 
force- free parameter a (where J = aB) which is 
constant along field lines. Figure [6] shows contours 
of a* found in this way for the initial state (top 
panel) and the final state (t = 290, lower panel). In 
both cases the sum of a* over all field lines in the 
domain is zero. In the final state the small scales in 
a* have been smoothed out leaving two large-scale 
patches of opposite signs. Note in addition that the 
extremes of the value of a* have been reduced. 

The two patches of a* shown in Figure [6] are sig- 
natures of the separation of the braided field during 
the relaxation into a simpler structure consisting of 
two weakly twisted magnetic flux tubes of oppo- 
sitely signed twist, embedded in an approximately 
uniform field. Here we provide two simple visualiza- 
tions of the field structure in order to demonstrate 
this process. In Figure [7] we show sample magnetic 
field lines in the initial (left) and final (right) states. 
In each case two sets of field lines are shown. First, 
a set of field lines is traced from three circles of ra- 
dius 0.1, 0.3 and 0.5 lying in the z = —24 plane and 
centered at x x = (1.415,-0.4875,-24). A second 
set of field lines is traced from three circles of the 
same radii centered at X2 = ( — 1.06,0.0375,-24). 
The centers xi and X2 of the two flux tubes have 
been chosen to lie at the maximum and minimum 
values of a* (at z = —24) in the final state. 

In the lower panels of Figure [7] a two dimensional 
representation of this separation is shown. We trace 
field lines from a regular grid of points (over the re- 
gion [x,y] E [±6, ±6]) on the lower boundary. In 
the figure we plot the intersections of these field 
lines with the upper boundary z = 24, colour-coded 
such that those field lines that begin at x > 
at z = — 24 are coloured white, while those with 
x < at z = — 24 are coloured black. The left 
panel demonstrates the complexity of the footpoint 
mapping generated by the braided field in the ini- 
tial state. The right panel shows the much simpler 
mapping obtained in the final state. 

In summary, we have found that in the final state 
of our simulation the magnetic field has two impor- 
tant properties: 



1. The magnetic field has been 'unbraided', be- 
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coming much more topologically simple. 



2. The field relaxes towards a non-linear, rather 
than a linear force- free field. 







Figure 7: Above: field lines traced from fixed lo- 
cations on the lower boundary (and coloured ac- 
cording to location on the lower boundary). Below: 
locations of intersection with the plane z = +24 of 
field lines traced from a regular grid in the plane 
z = —24. Field lines traced from locations x > on 
the lower boundary are coloured white, from x < 
black. Plots are made for t = (left) and t = 290 
(right) for the run with 77 = 10 -3 . 



4.2 Discussion 

The nature of the relaxed state following the resis- 
tive relaxation of astrophysical magnetic fields (such 
as found in the Solar corona) is an important and 
unresolved issue. A linear force-free field is often 
conjectured to be the end-state of general turbulent 
relaxation processes. This notion was put forward 



by Taylor (1974) to explain the turbulent relaxation 



of a reversed field pinch (RFP), and was applied to 



the corona by Heyvaerts and Priest (1984, 1992) 



Noting that the full set of invariants of ideal MHD 
do not hold in a resistive situation, Taylor suggested 
that as the magnetic field of an RFP decays toward 
a state of minimum energy, a single invariant - the 
total magnetic helicity (H = J A • B dV) - will 
serve as a constraint on the minimization. Formu- 
lating this as a variational problem the end state 
of relaxation is found to be a linear force-free field, 
which matches experimental data of the RFP rather 



well (e.g. Taylor, 1986). In the solar coronal case, 
on the other hand, it is not yet known whether the 
total helicity is the only constraint on relaxation. 
It is certainly true that there are fundamental dif- 
ferences between the the two plasma environments, 
the most obvious being that field lines in the solar 
corona do not close there, but rather can be consid- 
ered to be anchored in the solar photosphere. 
In a series of numerical simulations Browning and 



coworkers ( 


Browning and Van der Linden, 2003; 


[Browning et al. 


, 2008; 


Hood et al., 


2009) have fol- 



lowed the development of the kink instability in 
an initially non-linear force-free cylindrical loop. 
In their simulations, various initial distributions of 
a take the form of cylindrically symmetric step- 
functions with different profiles. The authors find 
the end-state of relaxation to be a linear force-free 
field. By contrast, Amari and Luciani (2000) - who 



also tracked a turbulent resistive evolution subse- 
quent to a kink instability - found as an end state 
a non-linear force-free field containing two distinct 
flux tubes. 

Hood et al. (2009|) have suggested that a smooth- 



ing of a over the full domain is natural in the pres- 
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ence of fully developed turbulence, i.e. with recon- 
nection regions sufficiently spatially overlapping (in 
the case in which a dominant vertical magnetic field 
B z is present this implies an overlapping of cur- 
rent sheets when projected along field lines onto 
the xy plane). However, as we proceed to demon- 
strate below, our study indeed satisfies such a cri- 
terion, though we do not obtain a linear force-free 
field, as discussed above. We have recently discov- 
ered one additional constraint on the resistive relax- 
ation of line-tied magnetic flux tubes, which demon- 
strates why a Taylor relaxation to the linear force- 
free state is prohibited in the braided field studied 
here (lYeates et al.l 120101). It is clear that further 



studies are necessary to determine the end state of 
relaxation even in isolated flux tubes or braids. 

To qualify the above discussion, it is worth em- 
phasising at this point that the final state of our 
simulations is not an exact force-free field, due to 
the finite plasma pressure and residual oscillations, 
as discussed above. Within each twisted flux tube, 
there is a balance between magnetic tension asso- 
ciated with their internal twist, and a combination 
of magnetic and plasma pressure. That the plasma 
pressure is enhanced in the central part of the do- 
main is not a surprise, since the plasma has been 
heated there by a multitude of reconnection events. 
Although the Taylor relaxation scenario assumes a 
very small plasma- /3, and we have used a plasma- f3 
of order 0.1, we do not expect that the plasma pres- 
sure plays a crucial role in determining the qualita- 
tive properties of the final state. We have repeated 
the simulation run with 77 = 10 -3 with the plasma- 
f3 reduced by a factor of 10 with no change in the 
qualitative results. 

5 Magnetic reconnection 
5.1 Introduction 

In order for the magnetic field to resistively relax, 
it is necessary that magnetic reconnection occurs, 
since a subset of the field lines in the volume have 
a non-trivial topology - see Figure [ija). Recon- 
nection allows a change of the magnetic topology, 
and thus release of the energy stored in the mag- 
netic field. It is becoming apparent that magnetic 
reconnection processes are fundamentally different 
in two dimensions (2D) and three dimensions (3D, 



clearly the regime applicable to the solar corona). 
In 2D, the rate at which the magnetic flux is trans- 
ferred from one flux domain to another (through a 
magnetic X-point) is given by the value of the elec- 
tric field evaluated at the null point. However, in 
3D, the electric field E and magnetic field B are no 
longer perpendicular, and it is the component of E 
parallel to B (denoted E\\) that is the crucial quan- 
tity. More specifically, magnetic reconnection can 
be defined to occur within a region in which 



-/ 



E\\ dsytO 



(7) 



where s is a parameter that runs along magnetic 
field lines, and where the quantity En should be 



spatially localised (Schindler et al. , 1988; Hesse and 



|Schindler , 1988). The rate of reconnection in this 
case is given by the maximum value of F over all 
field lines threading the diffusion region. Thus, the 
efficiency of the reconnection process is determined 
not by a local quantity evaluated at a given point, 
but rather by a 'global' property of the field in the 
vicinity of the diffusion region. Furthermore, the 
interpretation of this reconnection rate is far more 
complicated in 3D. The way in which the mag- 
netic flux evolves during the reconnection process 
has been shown to be fundamentally different in 2D 



and 3D (Priest et al. 


2003; 


Hornig and Priest , 2003 


Pontin et al., 2005a) 


. One crucial difference is that 



the magnetic field lines change their connectivity 
continuously throughout the diffusion region. A key 
result of this property is that magnetic field lines do 
not reconnect in a one-to-one 'cut and paste' fash- 
ion, as in simple 2D models. This has a profound 
effect on the way that the magnetic flux is restruc- 
tured by the reconnection process. The issue is fur- 
ther complicated by the fact that many different 
modes of magnetic reconnection may occur in 3D, 
depending on the local magnetic field structure in 
the vicinity of the diffusion region, with each hav- 
ing different characteristic behaviour in terms of the 



evolution of the flux (e.g. 


Hornig and Priest, 2003 


Pontin et al., 2004, 2005b; Wilmot-Smith and De 


Moortel , 


2007; Parnell et al. 


, 2010). 



5.2 Method of measurement 

Thus far we have described the formation of multi- 
ple current sheets during our simulations. We now 
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Figure 8: Plot of <E>(xo,?/o) (as defined in Equation 
^) for t = 50 for the run with 77 = 10 -3 . 

examine the properties of the magnetic reconnec- 
tion that occurs in these current sheets, in particu- 
lar the rate of the reconnection. Due to the strong 
background ('guide') field in the z-direction, there 
are no null points or closed field lines in our vol- 
ume. The magnetic reconnection takes place in the 
absence of such topological features, as described 



for example in the models of Hesse and Schindler 


(1988 


); 
); 


Priest and Forbes ( 


1992 


);] 


iornig and Priest 


(2003 


Wilmot-Smith et al. ( 


200( 


5, 2009c|). As dis- 



cussed above, the rate of reconnection is deter- 
mined, for an isolated diffusion region in 3D, by 
the maximum value of F = J E^ds attained along 
any field line threading the diffusion region. 



It is clear from Section [3.11 that there are multi- 
ple diffusion regions, at which reconnection is oc- 
curring, during the simulations. Therefore the first 
step towards calculating the reconnection rate is to 
first identify all of the individual reconnection re- 
gions that exist at each time. The global reconnec- 
tion rate is then obtained by summing the moduli 
of the maxima of F (see Eq.[7| over all such regions. 

At periodic intervals during the simulations, the 
data of all physical fields is output, providing a 
'snapshot' of the evolution. We proceed as follows, 
treating each individual snapshot in turn. 

1. Integrate E\\ along a set of field lines which in- 
tersect the z — plane in a rectangular grid 
with resolution approximately twice that of the 
computational mesh (to ensure that no infor- 
mation is lost in areas where the field line map- 
ping shows expansion). For clarity we inte- 
grate field lines only over the positive half- 



space, with the domain having symmetry about 
the z — plane. This yields a two-dimensional 
quantity 

s(z=2A) 

$(x ,yo)= J E\\ds (8) 

so 

subject to {x,y,z) — (xo,yo 5 0) at s — so, i-e. 
(#o,yo) i s the intersection of the appropriate 
field line with the z — plane. 

2. Search for the maxima and minima of <E>. We 
want to count only the reconnection rate for 
'isolated' diffusion regions, and so we require 
that $ falls below some threshold between 
peaks in order to class two regions as corre- 
sponding to separate reconnection processes. 
The threshold we choose to class two given 
maxima as 'separate' is that $ must drop be- 
low 60% of the lower of the two peak values 
between the peaks. (In principle we would like 
$ to fall to zero, but the effect of projecting 
along the field lines is likely to mean that re- 
gions that may be isolated in 3D overlap to 
some extent in the 2D projection. In addition, 
numerical quantities are never exactly zero in 
practice, so it is not realistic to enforce such a 
condition.) 

3. To be classed as a reconnection region, max- 
ima/minima are required to contain at least 
500 field lines. This corresponds to a cross- 
sectional area of around 0.02, or around 40 
grid cells. (This does not dramatically alter 
the reconnection rate obtained. However, it 
discounts large numbers of false positives in 
the algorithm used to identify isolated maxima, 
which otherwise selects many peaks which are 
sufficiently small (both in spatial coverage and 
modulus) to be close to the noise threshold.) 

4. To obtain the global reconnection rate, sum the 
moduli of the identified peak values, and multi- 
ply by 2 to obtain a value for the whole domain 
(since the above procedure is only implemented 
on the positive half-space in z). 

Below we present the results of implementing this 
procedure for the standard run with 77 = 10 -3 , be- 
fore going on in the next section to describe the ef- 
fect of varying the value of 77. As a frame of reference 
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for the numerical values obtained for the reconnec- 
tion rate, we note the following characteristic values 
in our simulations. Since p w 1 and the equations 
have been non-dimensionalised by setting fio = 1, 
we have |v^| « |B| « 1. The 'guide field' £> z has 
magnitude of order 1 for all time, and at t = for 
example we have peak values of B X1 B y of approxi- 
mately 0.3, and rms values of approximately 0.05. 

5.3 Nature and rate of reconnection pro- 
cess 

A sample $ plot is presented in Figure [SJ where the 
shading intensity shows the value of <3>. It is clear 
that there exist multiple distinct reconnection re- 
gions, which is consistent with the 3D visualisation 
of the current density presented in Figure [2| We 
now consider the qualitative evolution of $ in time, 
shown in Figure [9j For clarity only one sign of $ is 
shown (<E> < 0). At t — only a single reconnection 
region exists, which we see split and then intensify 
(t = 15) as time progresses, consistent with the for- 
mation of the initial pair of current sheets. As the 
evolution continues we see the formation of multiple 
distinct reconnection regions which fill the majority 
of the central portion of the xo2/crpl&ne where the 
twisted and braided flux is located (this filling is 
even more apparent when $ of both signs is plot- 
ted, see Figure[8]). As the simulation continues, and 
the current again assumes a large-scale structure as 
described above, the number of distinct reconnec- 
tion regions diminishes once more. 

Quantitative measures of the behaviour described 



above are presented in Figure [10| The global recon- 
nection rate - plotted as the solid curve in the left 
panel of Figure 10 - is defined by 2 x \&max\, as 
discussed above, where we multiply by 2 to take into 
account the lower half-space z < 0. The plot shows 
that initially the reconnection rate grows steadily, 
peaking at t = 35. There is then a gradual decay 
until t = 170 when the rate settles to an approxi- 
mately constant value (representing large-scale dif- 
fusion). The first thing to note is that the time of 
the peak reconnection rate (t = 35) does not coin- 
cide with the time of the peak current in the domain 
(t = 15). Furthermore, it does not occur at the time 
when the reconnection rate for any given single re- 
gion is a maximum (t = 15). Rather, the reconnec- 
tion rate continues to increase after this time, as a 
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Figure 9: Plots of $ for $ < only, at times t 
0, 15, 50, 140, 290, for the run with 77 = 10" 3 . 
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Figure 10: Left: Reconnection rate versus time 
(solid line) for the run with r\ = 10 -3 , defined 
by 2 x J2\<&max\- Also plotted (dashed line) is 
2 x (absolute maximum of |$|). Right: Cumulative 
reconnected flux versus time, for the same run. 



result of the fragmentation of the spatial location of 
the reconnection process. So the reconnection rate 
grows not because of an intensification of the cur- 
rent or an increase of the rate in any one region, but 
because of an increase in the fragmentation of the 
volume in which reconnection processes take place. 
This is demonstrated by the dashed line in the left 



panel of Figure [TUJ which shows the maximum rate 
for the fastest reconnection process (of either sign) , 
defined by \$\ global max, multiplied again by 2 for 
consistency. The number of individual reconnec- 
tion regions identified throughout the simulation is 
shown by the dashed line in the top panel of Figure 

The reconnection rate as calculated above is of 
course only an estimation. We would argue that it 
provides a lower bound on the reconnection rate for 
the system. There are a number of reasons why it 
should give a conservative estimate: 

• finite spatial resolution of the simulation grid, 
which will in general underestimate the true 
current maxima, 

• finite spatial resolution of the field line grid 
means the actual peak $ will not be exactly 
obtained, 

• finite temporal resolution will mean that we fail 
to capture the highest global peak value, 



along some field lines will degrade some max- 
ima, 

• exclusion of small regions (step 3 above), 

• projection along field lines may lead to regions 
distinct in 3D failing the test for separate re- 
gions, criterion described in step 2 above. 

5.4 Reconnected flux 

It is straightforward to estimate the total recon- 
nected magnetic flux from the reconnection rate 
obtained as described in the previous section, by 
performing a simple time integral. The cumulative 
reconnected flux is plotted against time in the right 
panel of Figure [TOj It turns out that the total num- 
ber of (non-dimensional) units of flux reconnected 
during the simulation is « 41.2. It is revealing to 
compare this figure with the total poloidal flux in 
the initial state, which is 30 units. (During the ideal 
relaxation used to define the initial condition for the 
present resistive MHD simulations, this quantity is 
conserved. In the symmetric pre-initial field, which 
is known in closed form, the poloidal flux is equiva- 
lent to the total flux of either sign passing through 
the y = plane.) Furthermore, as we saw in Section 
[1J the final state of the relaxation is not the homo- 
geneous field, but rather contains a finite amount 
of twist and thus poloidal flux. Due to the high 
degree of symmetry in the final state, the remain- 
ing poloidal flux may be estimated by integrating 
the positive (say) flux B through the y = plane. 
The result that we obtain is 15.3 units, indicating 
that around half of the initial poloidal flux has been 
cancelled. 

To recap then, beginning with 30 units of poloidal 
flux, 41.2 units are reconnected and 15.3 remain. 
So the total reconnected flux is greater than the 
quantity that would appear to be available. The 
only conclusion is that (at least some of) this flux 
is reconnected multiple times, most likely in differ- 
ent reconnection regions, during the relaxation and 
'unbraiding' of the field. On average, each unit of 
flux is reconnected around 2.8 times. This demon- 
strates the highly complex nature of the reconnec- 
tion process that occurs in complex (realistic) mag- 
netic fields such is the one considered here. 



The result should be compared with that of Par 



partial cancellation of positive/negative E\\ nell et al. (2008), who observed what they termed 
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'recursive reconnection' in a numerical simulation 
in which two opposite-polarity flux sources were 
driven past one another in the presence of an over- 
lying field. In their simulation, a number of sepa- 
rator lines (field lines connecting one magnetic null 
point to another) were formed through which the 
reconnection took place. Thus the reconnection re- 
gions were long-lived and individually identifiable 
from one time-frame of their simulation to another. 
They calculated that the magnetic flux was recon- 
nected on average 1.8 times, and argued that this 
occurred cyclically around the same circuit of re- 
connection regions (separators). It is clear that we 
observe a very similar phenomenon occurring here, 
even though the magnetic field configuration is com- 
pletely different (and, additionally, isn't continually 
being driven). The difference in the present case is 
that due to the lack of topological features present 
in the magnetic field and the high degree of com- 
plexity, it is difficult to label individual reconnection 
processes and follow them in time. The magnetic 
flux undergoing the reconnection process described 
above then may be best described as being 'mul- 
tiply' rather than 'recursively' reconnected. Given 
that such similar processes occur in such vastly dif- 
ferent configurations, this suggests that such 'mul- 
tiple' or 'recursive' reconnection is likely to be ubiq- 
uitous in astrophysical magnetic fields. 



6 Dependence on resistivity 

So far we have analysed the resistive relaxation 
properties of our model braided magnetic field for 
a fixed value of the resistivity. However, the dissi- 
pative length scale in the run we have so far anal- 
ysed is some orders of magnitude larger than length 
scales present in the complex field line mapping at 
t = (see |Wilmot-Smith et al.[ |2009a[). 



We there- 
fore proceed in this section to describe the depen- 
dence of our relaxation process on the value of the 
constant resistivity, 77. We have repeated our simu- 
lation at one higher and one lower value of 77, namely 
77 = 10 -2 and 77 = 2 x 10 -4 (the latter being the low- 
est value we can use at the present grid resolution 
while maintaining a reasonable resolution of current 
sheets in the system). 

Qualitatively, the differences between the simu- 
lation runs can be described as follows. The ini- 




Figure 11: Plots of $ for $ < only, at 
50, 140, 290, for the run with 77 = 2 x 10" 
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Figure 12: Number of identified reconnection re- 
gions (upper panel), reconnection rate (middle 
panel) and cumulative reconnected flux (lower 
panel) for runs with 77 = 10 -2 (grey solid line), 
77 = 1CT 3 (black, dashed) and 77 = 2 x 10 -4 (black, 
solid). 



tial magnetic field loses stability at a time that 
is independent of 77, with two current layers form- 
ing at locations that are also independent of 77 (see 
Paper I). However, the subsequent evolution does 
vary greatly between the simulation runs. First, 
the lower the value of 77 used, the greater the de- 
gree of fragmentation of the current structure. This 
is quantified below where we also discuss the re- 
connection rate and evolution of the magnetic flux. 
It can also be seen by comparison of the frames 
showing one sign of $ in Figures [9] and 11 The 
second main effect of reducing 77 is to increase the 
timescale of the relaxation process. In other words, 
the plasma resistivity acts to inhibit the tendency 
to form an increasingly complex system of current 
layers, such that for larger 77 the system transitions 
to a diffusive-timescale evolution in a shorter period 



of time (again, compare Figures [9] and 11). 

The above description can be quantified by 
analysing the global reconnection rate in the same 
way as before. First, in Figure [l2] (top panel) we 
plot the time evolution of the number of individual 
localised reconnection processes within the domain, 
for the runs with different 77. It is clear that for 
lower 77 not only is the peak number higher, but the 
period over which the reconnection process occurs 
in a fragmented volume is significantly longer. Of 
course our algorithm is not sufficiently robust for 
the exact numbers of regions identified to have any 
physical meaning - however, we are confident that 
the trends that are demonstrated when varying 77 
are physical. Next, in the middle panel of Figure 
[12] we compare the measured global reconnection 
rate versus time. The plots indicate that there is 
a dependence of reconnection rate on resistivity - 
albeit a rather weak one. However, once again the 
robustness of our algorithm for calculating <E> could 
be questioned here. It is also possible that if some 
small amount of numerical dissipation is present for 
the low 7] run, then we may under-estimate the re- 
connection rate by under-estimating 77. 

The greatest insight into the changing nature of 
the relaxation process is gained by plotting the cu- 
mulative reconnected flux versus time for the dif- 
ferent simulation runs - see the lower panel of Fig- 
Here we see that the lower 77 becomes, the 



ure 
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greater the total amount of flux reconnected in the 
cascade phase of our simulations. This in itself is 
not an absolute quantity (there is no precise cri- 
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terion that we can apply for the transition from a 
fast, 'reconnective' evolution to a slow diffusive evo- 
lution), but we can compare this number both with 
the initial poloidal flux and with the net twist re- 
maining in the field at the end of our runs. Note 
that in the final state of each simulation the field 
has the same structure as described in Section H] of 
two oppositely twisted flux tubes. In the run with 
77 = 10 -2 our measured cumulative reconnected flux 
is 24.5 units while 13.1 units remain (compare this 
with the 30 units of poloidal flux in the initial con- 
dition). As discussed in Section 5.4, for 77 = 10 -3 



we measure 41.2 units reconnected and 15.3 units 
remaining. Finally, for the run with 77 = 2 x 10 -4 
we measure 61.8 total units reconnected, while 16.2 
units remain. We observe that there is a clear trend 
as 77 is reduced for more reconnection to be required 
to 'unbraid' the field: the flux is reconnected on av- 
erage 1.4 times for 77 = 10 -2 , 2.8 times for 77 = 10 -3 
and 4.4 times for 77 = 2 x 10 -4 . We would expect 
this trend to continue, at least until 77 has a suf- 
ficiently low value that structures on the scale of 
the initial field line mapping are resolved. As noted 



by Wilmot-Smith et al. (2009b), these length scales 



decrease exponentially with increasing braid com- 
plexity. Therefore in a true braided coronal loop 
we would expect a highly complex network of cur- 
rent sheets threading the loop, in which the flux is 
multiply reconnected to a high degree. 



7 Conclusions 

We have described a series of numerical experi- 
ments in which we considered the resistive MHD 
evolution of a braided magnetic field between two 
perfectly conducting parallel plates. Although the 
magnetic field in the simulations is initially approx- 
imately force-free, it was described in Paper I that 
the field experiences an instability on some char- 
acteristic time-scale. Here we considered the sub- 
sequent evolution of the system, which is best de- 
scribed as a resistive relaxation. The route taken 
to find a new equilibrium involves the formation of 
a complex array of current sheets which are scat- 
tered throughout the domain. The current sheets 
have a ribbon-like appearance, tending to be highly 
elongated in the direction along the loop, i.e. paral- 
lel to the strong axial field. The formation of this 



myriad of current layers suggests that a turbulent 
cascade may develop for higher Reynolds numbers 
than we have been able to use. Due to the absence 
of boundary forcing, this turbulence would fall un- 
der the heading of decaying turbulence. It would 
be of interest to push the spatial resolution of our 
simulations (thus allowing the reduction of 77) to 
investigate whether a regime of fully-developed tur- 
bulence arises, by examining spectral properties of 
the magnetic and velocity fields. This is beyond 
the scope of this paper, and we leave it to a future 
investigation. 

The relaxation process as described above results 
in a simplification of the magnetic field structure 
as demonstrated by the mapping of field lines from 
one line-tied boundary to the other (Figure [7]). The 
magnetic field lines are untangled, such that no 
three field lines are braided about one another any 
longer in the final state. The final state in fact con- 
sists of two (weakly) twisted flux tubes, of oppo- 
sitely signed twist, embedded in an approximately 
uniform field. This final state approximates a force- 
free field, V x B = ckB, in which field lines in one 
twisted flux tube have positive a and in the other 
negative a. Thus the system approaches a non- 
linear force- free field, which is not consistent with 
the Taylor relaxation picture put forward by |Hey- 
vaerts and Priest (1984). In principle since the net 
current in the system is zero, a Taylor relaxation 
would lead to a final state with a = 0on every field 
line, i.e. the homogeneous potential field. Clearly 
the single constraint (global helicity conservation) 
of the Taylor hypothesis is not sufficient to describe 
this relaxation. Extra constraints on the relaxation 
which prevent some of the helicity cancellation must 
be present. Although it is always dangerous to place 
too much trust in the extrapolation of simulation re- 
sults to the extremely high Reynolds number coro- 
nal plasma, it appears unlikely that the cause could 
be related to the 'turbulence' in the system not be- 
ing sufficiently developed. In fact we retain more 
twist in our flux tubes in the final state for lower 77, 
where the turbulence is better developed. In fact, 
an additional constraint on the relaxation has re- 



cently been discovered (Yeates et al. , 2010). 



The 'unbraiding' of the magnetic flux during the 
relaxation to a non-linear force-free field occurs by 
magnetic reconnection. This reconnection occurs 
in the absence of magnetic nulls or closed field 
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lines (note the presence of a strong background 
field throughout the domain). Furthermore, recon- 
nection occurs in a multitude of regions that are 
spread throughout the volume. In order to deter- 
mine the efficiency of the reconnection process, we 
must therefore sum the reconnection rate over all 
reconnection diffusion regions within the volume. 
The outcomes of performing such an analysis are as 
follows: 

1. The global reconnection rate continues to grow 
for some time after the peak current in the do- 
main begins to fall. That is, during the inter- 
mediate stages of the simulations, the recon- 
nection rate grows not because of an intensifi- 
cation of the current or an increase of the rate 
in any one region, but rather because of an in- 
crease in the fragmentation of the volume in 
which the reconnection processes take place. 

2. The peak value of the global reconnection rate 
is at most weakly dependent on the resistivity. 

3. The number of identifiable reconnection re- 
gions increases as the resistivity is decreased. 

4. The total quantity of magnetic flux that is re- 
connected is greater than the total poloidal flux 
present. This implies that the magnetic flux is 
'multiply-reconnected' in the complex array of 
reconnection regions. 

5. As the resistivity is decreased - resulting in an 
increase in the complexity of the multitude of 
current layers - the average number of recon- 
nect ions for each unit of flux increases. 

In summary, the loss of stability and subsequent 
relaxation of braided coronal loops results in a com- 
plex array of current layers which permit a dissipa- 
tion of the stored magnetic energy and the attain- 
ment of a lower energy non-linear force- free field. 
These results yield a new and intriguing potential 
resolution to the long-standing debate over the va- 
lidity of Parker's model of coronal heating by topo- 
logical dissipation. Specifically, while current singu- 
larities are not an inevitable result of random field 
line tangling, once the field line mapping becomes 
sufficiently complex, an instability will be triggered, 
causing the electric current structures to collapse to 
small scales. The complexity that is inherent in the 



field then results in the formation of a fragmented 
system of current layers (the sum the associated re- 
connection processes being 'fast' in the sense that 
the rate of energy release is only weakly dependent 
on the magnetic Reynolds number). The global re- 
sult is a turbulent dissipation of the excess mag- 
netic energy stored in the field, ultimately in the 
form of heat. There are many interesting aspects of 
the above study that warrant further investigation, 
such as determining what parameters and physical 
quantities govern the nature of the final state, the 
possible turbulence properties of the relaxation, and 
the spatial and temporal distribution of the result- 
ing heating in the loop. 
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